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Abstract 

~ An  efficient  numerical  (finite  element)  method  is  presented  for  the 
dynamic  analysis  of  rapidly  propagating  cracks  in  finite  bodies,  of  arbitrary 
shape,  wherein  linear-elastic  material  behavior  and  two-dimensional  conditions 
prevail.  Procedures  to  embed  analytical  asymptotic  solutions  for  singularities 
in  stresses/strains  near  the  propagating  crack-tip,  to  account  for  the  spatial 
movement  of  these  singularities  along  with  the  crack-tip,  and  to  directly 
compute  the  dynamic  stress-intensity  factor,  are  presented-JJumerical  solutions 
of  several  problems  and  pertinent  discussions  are  presented  ita  Part  II  of 
this  paper.  ' 
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Introduction 


A  concise  summary  of  the  present  status  of  the  theories  of  dynamic  crack 
propagation  can  be  found  in  a  recent  article  by  Freund  [1].  Several  analytical 
solutions  of  the  linear  elasto-dynamic  equations  for  crack  propagation  in 
unbounded  plane  bodies  have  appeared  earlier.  These  include  the  works  of: 

Yoff&;  Cragg;  Broberg;  and  Baker,  for  Mode  I  (plane-strain  opening  mode) 
crack  propagation;  and  the  works  of:  Eshelby;  and  Achenbach,  for  Mode  III  crack 
extension.  All  the  above  works  are  summarized  and  referenced  in  a  paper  by 
Freund  [2],  who  considered  the  problem  of  a  hlaf-plane  crack,  in  an  elastic 
solid  subject  to  time-independent  loading,  which  is  initially  at  rest  and, 
at  a  certain  instant,  begins  to  move  with  either  a  constant  velocity  [2]  or  a 
non-uniform  velocity  [3].  The  studies  in  [2,3]  were  later  extended  [4]  to 
consider  stress-wave  loading.  However,  as  is  usually  the  case,  to  study 
dynamic  crack  propagation  in  finite  bodies  of  arbitrary  geometry,  it  is 
necessary  to  formulate  consistent  numerical  methods,  which  may  capitalize 
on  the  insights,  into  the  field  behaviour  near  propagating  crack-tips,  gained 
through  the  analytical  solutions.  A  critical  appraisal  of  several  and  varied 
numerical  solution  techniques  in  dynamic  fracture  mechanics  was  made  in  a  1978 
paper  by  Kanninen  [5].  Most  of  the  dynamic  finite  element  methods,  for  fast 
crack-propagation  analysis,  reviewed  in  [5J  use  the  conventional  finite  elements 
with  simple  polynomials  for  assumed  displacements,  and  do  not  account  for  the 
singularity  in  strains  near  the  crack-tip.  Further,  in  these  methods,  the 
dynamic  crack  propagation  was  simulated  by  a  "gradual"  release  of  the 
restraining  nodal  force  at  a  finite  element  node  which  represents  the  "current" 
crack-tip.  The  dynamic  stress-intensity  factor  is  then  extracted  from  the 
displacement  field  or  from  the  work  done  in  releasing  the  nodal  force.  It 
was  concluded  in  [5]  that  the  above  "node-release"  techniques  were  not  sufficient 
ly  accurate. 
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Since  the  appearance  of  [5],  Bazant  et  al  [6]  have  presented  a  calibrated, 
non-singular,  crack-tip  element  procedure  for  the  dynamic  analysis  of  running 
cracks.  In  the  procedure  of  [6],  the  finite-element  grid  moves  undeformed 
with  the  crack-tip.  However,  the  procedure  of  [6]  has  two  serious  limitations: 
(i)  it  is  restricted  to  finite  bodies  whose  surfaces  and/or  bimaterial  inter¬ 
faces  are  parallel  to  the  direction  of  crack  propagation;  and  (ii)  more 
importantly,  it  cannot  be  applied  to  bodies  having  finite  dimensions  in  the 
direction  of  crack  propagation.  On  the  other  hand,  Aoki  et  al  [7]  presented 
a  finite  element  procedure  wherein  the  singular  nature  of  stress/strain 
near  the  propagating  crack-tip  is  accounted  for  a  priori.  However,  in  [7], 
only  when  the  crack-tip  has  reached  close  to  the  boundary  of  the  singular 
element,  the  entire  singular  element  is  shifted,  as  a  rigid  body,  to  a  new 
location.  The  numerical  details  of  the  procedures  are  still  somewhat  sketchy 
in  [7].  Finally,  King  and  Malluck  [8]  reported  a  procedure  of  simulating 
crack-propagation  similar  to  that  in  [7],  except  that  the  singular-element  used 
in  [8]  has,  built  w'thin  it,  a  large  number  of  eigen-function  solutions 
corresponding  to  a  stationary  crack.  In  an  attempted  simulation  of  the  well- 
known  problem  of  Baker,  the  procedure  in  [8]  produced  spurious  oscillations, 
of  large  amplitude,  in  the  solution  for  dynamic  stress-intensity  factor,  as 
compared  to  the  analytical  solution.  Based  on  these  results,  it  is  suggested 
in  [8]  that  the  procedure  in  [8]  may  not  be  feasible  for  simulating  large 
scale  fast  fracture. 

In  Part  I  of  the  present  paper,  a  "moving  singular-element"  procedure  is 
presented  for  the  dynamic  analysis  of  fast  crack-propagation  problems  in 
arbitrary  shaped  finite  bodies.  In  the  present  procedure  a  singular-element, 
within  which  a  large  number  of  analytical  eigen  functions  corresponding  to  a 
propagating  crack  are  used  as  basis  functions  for  displacements,  may  move  by  an 
arbitrary  amount  AL  in  each  time- increment  At  of  the  numerical  tine-integration 


3 


procedure  (as  opposed  several  time  steps,  say  6  to  8,  per  increment  of  crack 
growth,  used  in  the  procedures  reviewed  in  [5]).  The  moving  singular-element, 
within  which  the  crack-tip  always  has  a  fixed  location,  retains  its  shape  at 
all  times,  but  the  mesh  of  "regular"  (isoparametric)  finite  elements,  surround¬ 
ing  the  moving  singular-element,  deforms  accordingly.  An  energy-consistent 
variational  statement  is  first  developed,  as  a  basis  for  the  above  "moving 
singular-element"  finite-element  method  of  dynamic  crack  growth  analysis. 

The  present  procedure  leads  to  a  direct  evaluation  of  dynamic-stress  intensity 
factor(s),  since  they  are  unknown  parameters  in  the  assumed  basis  functions 
for  the  singular-element. 

In  Part  II  of  the  paper,  several  numerical  results  for  cracks  propagating 
in  finite  bodies  are  presented  and  discussed. 

In  the  following  we  discuss  the  details  of  formulation  of  a  moving-singularity 
finite  element  formulation  for  analyzing  dynamic  crack  propagation. 


I.  Basis  Functions  for  a  Moving  Singular-Element 

We  consider  Mode  I-type  dynamic  crack  propagation  in  two-dimensional 
(plane  strain)  linear  elastic  isotropic  bodies  of  finite  geometry.  Let  xa 
(a  *  1,2)  be  fixed  cartesian  coordinates  in  the  plane  of  the  body,  and  x^  be 
the  thickness  coordinate  of  the  body  such  that  X2  *  0  defines  the  plane  of 
the  crack.  In  the,  context  of  the  present  numerical  method,  without  loss  of 
generality,  we  consider  the  case  when  the  crack-tip  is  moving  along  x^  axis 
at  a  constant  speed  v.  We  introduce  the  coordinate  system  (£,  X2)  which  remains 
fixed  with  respect  to  the  moving  crack-tip,  such  that  £  *  x^  -  vt.  Let  $  and 
4*  be  the  dilatational  and  shear  wave  potentials,  respectively;  and  let 
and  Cg  be  the  corresponding  wave  speeds.  It  can  then  be  shown  [2]  that  is 
governed  by  the  equation: 
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(1.1) 


II  -  <v/c)2i  & + £* 
d  3C2  3x2 


-(2v/cJ) 


32»  +  a/cb  & 

3t35  3,2 


and  that  ¥  is  governed  by  a  similar  equation,  except  that  is  to  be  replaced 
by  C  .  Consider  the  "steady-state"  solution  to  the  homogeneous  part  of  the 
above  equation,  that  is,  the  solution  which  appears  time- invariant  to  an 
observer  moving  with  the  crack-tip.  This  eigen-function  solution  which  satisfies 
the  traction-free  condition  on  the  crack  face  m  +0) ,  can  be  derived 

easily,  as  for  instance  in  [9,10],  and  is  given  in  Appendix  A  for  the  sake 
of  completeness. 

In  the  present  procedure,  a  finite  region  (which,  for  convenience,  is 
taken  to  be  rectangular  in  shape)  near  the  moving  crack-  tip  is  modeled  by  one 
finite-element,  in  which  the  displacement  field  is  assumed  to  be  a  linear- 


where  (~)  and  (:)  under  a  symbol  denote  a  column  vector  and  a  matrix,  respec¬ 
tively;  and  uS  denotes  the  vector  of  displacements  in  the  singular-element. 

We  note  that  the  total  velocity  and  acceleration  of  a  material  point  in  the 
singular-element  are  given  by: 


and 


u8  «  U  8  -  v(U)  g 

*w  ***  *•  *•  c  ^ 

uS  -  U  B  -  2v(U),£  0  +  v2(U),^  0 


(1.5) 

(1.6) 


where,  a  (*)  denotes  a  total  derivative  with  respect  to  time  t,  and  (  ),^ 

denotes  a  partial  derivative  with  respect  to  £. 

Let  the  domain  of  the  singular  element  in  the  present  procedure  be  V 

s 

and  its  boundary  be  3V  ;  and  let  p  be  that  part  of  3V  where  the  usual  iso- 

parameteric  finite  elements  adjoin.  In  order  that  convergence  of  the  present 

finite  element  method  may  be  achieved,  compatibility  of  displacements,  velocities, 

and  accelerations  between  the  singular  elements  and  surrounding  regular  elements, 

i.e.,  at  Pg,  is  maintained  in  a  least  squares  sense  as  described  below.  Let 

the  displacement,  velocity  and  acceleration  assumption  for  the  regular  element, 

at  p  ,  be  taken,  respectively,  as: 
s 


u 


R 


(I.7a,b,c) 


where  N  are  functions  of  the  boundary  coordinate  n(xa)  at  3Vg,  and  qg  is  the 

•  •• 

vector  of  displacements  at  nodes  at  pg.  The  parameters  B,  B  and  B  are  so 
chosen  that  they  minimize  the  error  functionals : 


R.2 

u  )  dp 


.  R.  2  , 

-  u  )  dp 


Using  Eqs.  (I. A, 5, 6  and  7)  in  (1.8-10),  and  minimizing  1^,  It 

•  •• 

ly  with  respect  to  0,  3  and  B  it  can  be  shown  that. 


..R.  2  . 

-  u  )  dp 

(1.3.9,10) 
and  I.j  success  ive- 
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B  **  A  q  i  S  “  A  5  ^  B  q  >  3  “  A  q  ^  ^  B  q  ^  C  q  (X. 11, 12, 13) 

*“  •»  ~ S  «  <v  <o*S  o#  -vS  5j  —8  oj  48  •»  48 


where : 


=  H-1  G  ; 

B  *  (v)  H_1  E  A 

ti.i6;i5) 

-  2(v)  H-1  E 

B  -  (v2)  H"1  F  A 

«w  >w  •* 

*v 

(1.16) 

*  f  UT  U  dp 

•V  = 

;  G  -  f  UT  N  dp 

**  •'D  *  ~ 

HS 

(I.17a,b) 

-If  <?•* 

^s 

ip  ;  F  -y*  UT  (U),^dp 
^s~ 

(I.17c,d) 

Thus,  Eqs.  (1.4, 5, 6)  together  with  (1.11,12,13)  represent  the  displacements, 
total  velocities  and  total  accelerations  in  the  singular  element,  in  terms 
of  its  nodal  displacements,  velocities,  and  accelerations,  q  ,  q  and  q  , 
respectively.  Thus,  if  q  at  p  is  determined,  then  3  (and  especially  the 
mode-I  stress  intensity  factor  B^),  can  be  determined  directly.  Finally, 
it  is  noted  that  the  above  Eqs.  (1.4, 5, 6)  and  (1.11,12,13)  represent  the 
assumptions  for  the  relevant  field  variables  in  the  singular  element  at  any 
generic  time  t. 

Now  we  consider  the  problem  of  dynamic  crack  propagation  within  a  time 
increment  At  between  two  generic  times  t^  and  t2* 


I.B.  Variational 'Principle  for  Dynamic  Crack  Propagation  Analysis 

In  the  following,  we  present  a  variational  statement  for  dynamically 
growing  cracks  in  linear  elastic  solids.  Consider  two  instants  of  time  t^ 
and  tj  (=  t^  +  At)  at  which  the  variables  of  the  problem  are  denoted  by  super¬ 
scripts  1  and  2,  respectively.  At  time  t^,  let  the  volume  of  the  solid  be 
V^,  the  external  boundary  of  the  solid  where  tractions  are  prescribed, 

be  S^;  and  let  E*  and  E^  be,  respectively,  the  two  surfaces  of  the  crack. 

-2 

Also,  let  be  body  forces  per  unit  volume  In  the  body  at  time  t2» 


-  Am 
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We  assume  that  between  time  t^  and  t2»  the  crack  surfaces  change  by  AE. 


r 


o 


The  orientation  of  AE  to  E,  can  be  determined  by  some  crack-growth  direction 

criterion;  however,  for  pure  Mode  I,  self-similar  growth  is  assumed.  The 

newly  created  crack-surfaces  can  be  traction  free,  but,  for  the  sake  of 

*  -2+  -2- 

generality,  assume  that  new  tractions  T^  and  T^  are  applied  on  the  new  crac.: 

+  —  —2 
faces  AE  and  AE  ,  respectively;  likewise,  let  new  tractions  T^  act  at  So2. 

The  principle  of  virtual  work  applied  at  t2  can  be  written  as: 


0  +  p“i  dv  ~  J  ^i  ^ui  "  ^  ^i  5ui  ds 

2  V2  Sa2 

+  <?i)+  <S“i>+ds  -L  <*!>'  d‘ 

L1  \ 

- 

-I  .  (T?)+  (6u^ds  -  (T h~  (6u*)"  ds 

•'Ar  1  1  Ar“  1  1 


(1.18) 


However,  for  the  case  of  cracked  structures,  the  changes  in  volume  and  external 
surfaces  between  times  t^  and  t2,  due  to  a  change  in  the  crack-surface  by 
AE  alone,  can  be  assumed  to  be  negligible,  i.e.,  -  V2  and  So1  =  S^*  It 

2+  2—  2  +  2  - 

is  important  to  note  in  Eq.  (1.18)  that  (u^)  #  (u^)  [or  (<5u^)  r  (&k)  ]  at 

the  initial  crack  surfaces  E*  and  E^,  nor,  more  importantly,  for  the  newly 

created  crack  faces  AE+  and  AE  during  the  time  interval  t2  -  t^  (“At).  If 

2  2  +  2  — 

similar  virtual  displacements  (<5u^),  such  that  (<5u^)  /  (6u^)  either  on  E, 

or  on  AE,  are  considered  in  the  statement  of  virtual  work  at  time  t^  (prior 
to  the  creation  of  new  crack  faces  AE) ;  this  statement  can  be  written  as: 


ij 
■  .1 

1 


*It  is  noted  that  the  element  basis  functions  assumed  in  Eqs.  (1.2)  and  (1.3) 
satisfy  only  the  traction-free  conditions  on  the  crack-face.  It  is,  however, 
easy  to  accommodate  non-zero  traction  conditions  on  the  crack-face  by  introduc¬ 
ing  appropriate  additional  terms  in  Eqs.  (1.2  and  1.3).  These  additional 
terms  are  so  chosen  that  they  satisfy  the  non-zero  crack-face  traction  con¬ 
ditions  either  exactly  or  in  an  average  sense. 


■i 


•  .  ~r* 


1 
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(1.19) 


0  -  I  (o±J  <5e±j  +  pu±  ou^)  dv  -  /  F±  6u±  -  /  Tj  6u‘  ds 

2  v.  SJ 

°2 

'/+  (5i)+  (Su2)+ds-/*_  (fh"  (6u2)'ds  -  /*  .(aj,  vj)+  (Su2)+  ds 
1  1  •'Z  1  i  */AZ+  iJ  J  1 

•  /  <aL  vh"  (6u.2)"  ds 

wherein  the  approximations  V£  «  ~  are  used  and  is  a  unit  normal 

to  Z^.  Adding  Eqs.  (1.18)  and  (1.19),  the  virtual  work  principle  governing 
dynamic  crack  propagation  between  times  t^  and  t£  can  be  written  as: 

/  {(0ij  +  °ij)6eij  +  p(Ui  +  Ui>  6ui  r-  (?i  +  fi)6ui}dv 

V2 

•f  (il +  Suids  *L  (ii +  fi>+  <6“i>+  *> 


+/-  +  *i>‘  +/  +  «i  +  vj)+  (6»*>+ 


+  /  (T2  +  oj  vj)“  (6u?)“  ds 

AZ~  3 


(1.20) 


In  the  finite  element  development,  the  domain  can  be  considered  to  be  broken 
into:  a  singular  element  surrounding  the  crack  tip  (See  Fig.  1),  and  a 

number,  N,  of  regular  elements  (n=l..N)  (thus,  V2  ’  V2s  +  5  V2Ri?  ;  likewise 

'  4-  +  + 

S  ~  *  Z  S_nn  .  Also,  as  seen  from  Fig.  1,  Z-r-  **  Z-,  +  Z  Ztt_  .  Henceforth, 
o2  n  o2Rn  1  si  n  lRn 

for  simplicity,  we  use  symbols  Vs>  Vj^,  S^,  Z^,  and  Z— ^  instead  of  V2s’ 

V2Rn’  ^02Rn’  ^1’  and  ^TRn’  resPecC^-ve^y •  We  now  restrict  our  attention  to  the 
mode-I  case  only,  i.e.,  when  the  applied  loading  is  in  a  direction  normal  to 
the  crack  plane  and  is  symmetric  with  respect  to  the  crack  plane  for  all 


times  t.  Thus,  for  the  mode  I  case,  using  the  above  notation,  the  virtual  work 
equation  as  applicable  to  a  system  of  finite  elements  may  be  written  as: 


iff  +  p<“i +  =i>  -  <?!  ♦  ?i>  «»j}dv 

*-VRn 

*4  (,1  +  5“ldS  -/E+  <*i  +  *i>+  IE]  ♦  S  V  +  dv 

^  s 

p  (“i  +  “i*  6ui  dv}'~^.+  (*i  +  *i>+  6ui+  dl 


-4+  (i* +  a«  vj1)+  6u'+  ds " 0 


(I.2D 


Assuming  that  crack  growth  occurs  between  times  t^  and  t^  (which  can  be  deter¬ 
mined  by  an  appropriate  criterion,  in  the  so-called  "application"  calculations 
using  the  given  material  dynamic  fracture  toughness  as  an  input;  or  is  known, 
a  priori,  in  the  so-called  "generation  phase",  i.e.,  in  the  case  of  simulation 
of  known  crack-tip  time  history  data),  the  singular-element  is  translated, 
in  the  mode  I  case,  along  the  original  crack  axis,  by  an  appropriate  distance 
AE  from  its  location  at  time  t^,  as  shown  in  Fig.  1. 

It  is  important  to  note  that  in  the  present  procedure,  this  amount  AE 
is  not,  in  any  way,  related  to  the  distance  between  any  two  adjacent  finite 
element  nodes  at  time  t^;  as  is  the  case  with  most  common  finite  element  methods 
which  use  the  "node-release"  technique  in  the  simulation  of  dynamic  crack 
propagation.  As  can  be  seen  from  Fig.  1,  as  the  singular-element  is  translated 
by  AE  between  t^  and  t^,  the  nodal  pattern  of  the  surrounding  regular  elements 
also  changes  between  t^  and  t2*  It  is  to  this  readjusted  finite  element  mesh 
at  time  t 2  that  the  virtual  work  equation  in  Eq.  (1.21)  is  understood  to  be 
applied.  However,  it  is  also  noted  that  only  the  nodes  of  the  elements 
immediately  surrounding  the  singular-element  are  readjusted  due  to  crack- 
growth  of  amount  AE  between  t2  and  t^.  Thus,  one  has  to  obtain  data,  such 
as  displacements,  velocities,  and  accelerations,  at  time  t^,  at  the  new 
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nodes  of  the  regular  elements,  which  are  indicated  by  solid  circles  in  Fig.  1. 
This  data  can  be  determined,  using  elementary  interpolation  techniques,  from 
the  known  data,  at  time  t^,  at  the  ’old*  nodes  at  time  t^,  which  are  indicated 
by  open  circles  in  Fig.  1.  The  details  of  these  interpolation  techniques  are 
omitted  for  simplicity  and  will  be  reported  elsewhere.  Thus,  one  is  in  a  posi¬ 
tion  to  know  the  relevant  data  at  time  t^  at  new  nodes  and  (hence  new  elements) 
corresponding  to  the  mesh  in  t2 ;  and  to  assume  the  appropriate  basis  functions 


for  the  relevant  variables  at  time  t2  for  the  mesh  at  time  t2,  as  follows: 
Known  at  t^  for  the  mesh  at  t^: 


in 

V 

Rn’ 

2i = 

N  qx  ; 

£l  =  B  2i  ; 

Si  “  H  5  Si 

(1.22-24) 

“1 = 

N  qx  J 

si  -  ;  !i 

(1.25,26) 

in 

v 

Si  = 

!i!i  ; 

~i  *  Hi  h; 

-  viHi,s  !i 

(1.27,28) 

h = 

Si  !i- 

2  v  U.  c  + 

^  ^ 

V1  Hl,«  ~1- 

(1.29) 

Si' 

Ssi  5i 

;  °~i  *  ?i  5i 

»  Ii  =  ?1  §1 

(1.30-32) 

Assumed 

at  time  t 

t  for  the 

mesh  at  time 

in 

VRn: 

“2  = 

=  N  q2  ; 

52  *  5  32  ; 

?2  =  5  5  32 

(1.33-35) 

~2  : 

=  N  q2  ; 

32  *  ~  32 

(1.36,37) 

in 

V 

~2 

■  Hz  h 

•  Hz  "Hz  Hz 

-  VZ  Hz, 5  !z 

(1.38,40) 

S  2 

■  Hz  !z 

-  2  vzHz,?  Hz  ' 

»  ''zHz.ES  5z 

(1.41) 

52  ' 

■  h  5z 

•  Hz  ■  Hz  !z 

’  T2  ~  ?2  ^2 

(1.42-44) 

where  the  familiar  vector  representations  for  displacements,  strains,  stresses, 
and  tractions,  are  employed  as  u,  £,  o  and  T  respectively.  Also,  v^  and  v, 
are  velocities  of  the  crack-tip  at  times  t^  and  t2  respectively,  and  the  eigen 


functions  IK  and  U9  depend  on  v  and  v9  respectively. 


Using  Eqs.  (1.22-44)  in  Eq.  (1.21),  the  finite  element  equations,  for 
arbitrary  variations  5q2  and  6&2  can  be  written,  as  shown  in  Appendix  B,  as: 


K  32  +  ~  §2  *  §2  +  Si  "  £  Si  "  m  ^or  VRn  in  V2  “  Vs 


*  *  .  *  * 

K  q  _  +  D  q  ,  +  m  q~“Q  for  V 
~s  ~s2  ~s  ~s2  ~s  2s2  2s  s 


(1.45) 

(1.46) 


where  K,  m,  Q-,  Q  ,  K  ,  D  ,  and  m  are  defined  in  Appendix  B,  from  which  it 

~  2  UjS  %S  ^  S 

*  4c 

can  be  seen  that  the  metrices  K  and  D  are,  unfortunately,  unsymmetric, 

^  S  ^  S 

while  the  others  are  all  symmetric.  In  Eq.  (1.45),  q^  and  q^  are  displacements 
and  accelerations  at  t2  at  nodes  everywhere  in  and  at  the  boundary  of  the 
region  O^-vp ;  whereas,  qg2,  qg2,  and  qg2  are  displacements,  velocities, 
and  acceleration  at  t2  at  nodes  along  the  boundary  3vg  of  the  singular  element. 
When  Eqs.  (1.45,46)  are  assembled,  it  can  be  seen  that  the  resulting  global 
"stiffness"  and  "damping"  (which,  however,  is  not  a  physical  damping  term) 
matrices  have  only  a  "small"  degree  of  unsymmetry,  confined  to  those  rows  and 
columns  corresponding  to  nodes  around  the  singular  element.  We  can  use  the 
common  time-integration  schemes  to  integrate  Eqs.  (1.45-46).  In  particular, 
we  use  the  Newmark's  method  which  can  be  characterized  by  the  approximations: 


q2  =  (q2  -  qx)  -  C2  -  C-j  q^  (1.47) 

q2  -  (q2  -  qx)  -  C5  qx  -  ^  (1.48) 


where  =  (5/1  At)  ;  C2  =  (5/Y)  -  1  ;  C3  -  (^p)  [(5/Y)  -  2] 

C4  -  1/Y  (At)2  ;  C5  =»  1/ (Y At)  ;  =  Off )  -  1  (1.49) 

where,  in  Che  present  calculations,  Y  *  4,  5  =  ‘2  arc  used.  With  the  difference 

approximations  in  Eq.  (1,47,48),  and  similar  ones  for  a  ,,  and  q  ,  we  reduce 

~s2  ~s2 

Eqs.  (1.45,46)  to: 
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*»> 


«<twt  i. 


K  ,2  -  9  for  V2  -  Vs 


Js  5.2  '  3s  f°r  V. 


where 


K  -  K  +  C4  m 

Q  =  Q  +  Q.  -  K  qx  -  m  Ji  +  m  (C4  qx  +  Cj  51  +  Cft  §]) 


A  ^  £  £ 

K«.  -  K  +  C.  m  +  C-  D 
~®  ss  4  =s  1  ~s 


(1.50) 

(1.51) 

(1.52) 

(1.53) 

(1.54) 


Ss  '  Ss  +  ‘C4  5.1  +  cs  3,1  +  c6  3sl>  +  3.  <C1  3sl  +  C2  3,1  +  C3  3sl> 

(1.55) 

A  A 

where  K  is  symmetric;  however,  Kg  is  unsymmetric.  When  Eqs.  (1.50,51)  are 
assembled,  we  obtain,  the  final  algebraic  equations: 


IK*]  {qj}  -  {Q  } 


(1.56) 


where  the  stiffness  matrix  in  Eq.  (1.56)  is,  in  general,  unsymmetric,  but  the 
unsymmetry  is  confined  mainly  to  the  rows  and  columns  corresponding  to  nodes 
around  V  .  A  rather  simple  technique  of  iterative  solution  of  the  above 
equation,  based  on  the  decomposition  of  the  stiffness  matrix  into  symmetric 
and  skew-symmetric  parts,  as  below,  was  used. 

%  [K*  +  K*T]  {q*P}  =  {Q*}  -  Jg  [K*  -  K*x]  (q* (p—1) }  (1.57) 

for  any  pth- iteration.  In  all  the  solutions  obtained,  only  two  iterations 
were  found  to  adequate.  Once  q2  Is  computed  from  Eq.  (1.56),  the  solution  for 
time  t2  +  At  can  be  repeated,  with  the  approximations  for  the  initial  values 
q2  and  q2  as: 


ft 

32  ’  C4 

*  *  .  * 

*32  ~  3i*  ~  c5  3i  c6  3i 

(I. 58) 

.*  .* 
32  =  31 

+  c?  q*  +  C3  q* 

(1.59) 
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where  C^,  C^,  are  defined  earlier,  and  ■  At  (1-6);  and  Cg  *  6At  (where 

a  value  of  6  *  \  is  used  presently) . 

* 

Once  the  nodal  displacements  (and  hence  the  corresponding  displacements 

at  the  nodes  of  the  singular-element),  at  time  tj,  are  computed  from  Eq.  (1.57), 
the  unknown  parameters  $  (and  hence  the  dynamic  stress-intensity  factor  8^) 

in  the  singular-element  can  be  computed  from  Eq.  (I. 11). 

Using  Eqs.  (1,58  and  59)  as  initial  data,  the  time-integration  between 
the  time  steps  t2  and  t^  (t2  +  At)  can  be  carried  out  andj  thus,  the  process 
can  be  repeated  for  all  subsequent  time  intervals.  The  successive  growth  of 
the  crack,  for  a  representative  problem  is  schematically  illustrated  in  Fig. 

2. 

\ 

From  the  example  given  in  Fig.  2,  it  is  seen  that  the  singular-element 
(A)  remains  its  shape  at  all  times  but  the  regular-elements  (B)  in  the 
"immediate  surrounding"  of  the  singular-element  continually  distort.  How¬ 
ever,  in  the  above  example,  at  t  »  2.0  Usee,  elements  B  have  distorted 
sufficiently  so  that  the  use  of  isoparametric  approximations  in  these 
elements  may  introduce  spurious  numerical  errors.  For  this  reason,  as 
typified  by  the  above  example,  at  t  ■  2.0  us,  the  regular  elements  B  are  re¬ 
adjusted  as  shown  in  Fig.  2.  This  involves  a  simple  reinterpolation  of  data, 
in  'B'  type  elements  from  t  -  2.0  -  0  Us  to  t  -  2.0  +  0  us,  the  details  of 
which  are  omitted  for  brevity.  Finite  element  calculations  detailed  earlier 
can  be  repeated  for  the  readjusted  mesh  at  t  -  2.0  +  Ous  until  the  B  type  ele- 

T 

ments  become  so  distorted  that  another  readjustment  may  be  warranted.  These 
mesh  readjustments  were  found  to  be  easy  to  accomplish  in  the  computer  coding 
based  on  the  present  approach. 

Finally,  it  may  be  of  interest  to  note  that  in  the  present  singular 

it 

element,  19  eigen  functions  for  a  propagating  crack  (See  Appendix  A)  were  used 

*The  number  of  eigen  functions  plus  the  number  of  rigid  modes  must  be  greater 
than  or  equal  to  the  number  of  degrees  of  freedom  at  the  boundary'.  A  study 
of  the  effect  of  the  number  of  eigen  functions  used,  on  the  results  was  con¬ 
ducted,  by  varying  this  number  from  17  to  25.  The  results  varied  only  in¬ 
significantly  (ie.,  less  than  0.4X),  and  the  number  of  eigen  functions  was 
chosen  to  be  19  in  all  subsequent  computations. 


along  with  a  rigid-body  translation  mode  In  xx  direction;  whereas,  there  are 
18  degrees  of  freedom  along  the  boundary  pg  of  the  singular  element.  The 
regular  elements  were  of  the  common  8— noded  isoparametric  type. 

It  should  be  remarked  that  the  problems  dealt  with  in  the  present  paper 
are  limited  to  the  case  of  determining  the  stress-intensity  factor  at  the 
crack-tip  which  is  propagating  with  a  prescribed  velocity-time  history.  Thus 
the  presently  treated  problem  may  be  considered  to  fall  in  the  category  of 
"generation  phase  calculations"  in  the  -sense  defined  in  [5  ].  The  present 
procedure  may  be  used  to  simulate  the  experimentally  determined  crack-velocity¬ 
time  history  in  test  specimens,  such  as  the  double-cantilever-beam  (DCB) 
specimen  |11  ],  to  determine  the  velocity  dependent  dynamic  fracture  tough¬ 
ness.  Using  this  as  input  data,  the  problem  of  determining  the  crack-tip 
motion  in  plane  elastic-bodies  subject  to  Mode-1  type  dynamic  transient 
loading  may  be  treated.  This  second  phase  of  research,  which  is  the  so-called 
"application  phase"  in  the  sense  defined  in  15  ],  is  currently  being  completed, 
and  will  be  the  basis  of  a  forthcoming  paper. 

Finally,  we  wish  to  note  that  once  the  basic  features  of  the  procedure 
based  on  the  present  moving  singular-element,  with  embedded  propagating- 
crack  eigenfunctions,  are  well  understood,  the  numerical  procedure  can  be 
further  simplified.  This  can  be  accomplished,  for  instance,  by  using  the 
well-known  distorted  isoparametric  elements  (the  so-called  "quarter-point 
elements")  [12]  in  place  of  the  present  singular  element.  Eventhougli  the 
results  from  the  use  of  a  quarter-point  element  are  not  expected  to  be  as 
accurate  as  from  the  use  of  the  present  singular-element;  such  results,  witn 
a  suitable  calibration,  may  be  used  in  analyzing  large-scale  fast  fracture 
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in  practical  situations.  The  results  from  the  use  of  a  quarter-point 
element,  and  their  comparison  with  those  reported  in  Part  II  of  this  paper 
(using  the  present  singular-element),  will  be  reported  on  shortly.  Also, 
since  it  is  known  [10]  that  the  eigen-functions  for  a  crack  propagating  at 
constant  velocity  differ  significantly  in  their  behavior  from  those  for 
stationary  crack  only  at  very  high  speeds  (v  =  C  )  of  propagation,  the  present 
procedure  can  be  simplified,  for  practical  purposes,  by  using  the  stationary- 
crack  eigen-functions  in  the  singular-element.  The  results  from  this  modifica¬ 
tions,  are  also  to  be  reported  shortly. 

Closure 

In  this  paper  we  have  presented  a  new  translating-singularity  finite 
element  procedure,  wherein  use  is  made  of  analytical  eigen-functions  for  a 
two-dimensional  crack  whose  tip  propagates  at  a  constant  velocity.  The  pro¬ 
cedure  is  capable  of  modeling  large-scale  fast  crack  propagation  in  finite 
two-dimensional  bodies  of  arbitrary  shape.  However,  the  type  of  problems 
considered  is  limited  to  the  case  of  determing  the  dynamic  stress-intensity 
factor  at  the  crack-tip  which  is  propagating  with  a  prescribed  velocity-time 
history. 

Implementation  of  the  present  approach  and  numerical  example  are  discussed 
in  an  accompanying  Part  II  of  the  paper. 
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Details  of  basis  functions  for  the  singular  element ,  for  the  mode  I  case, 
are  given  here.  The  eigen-functions  given  here  are  solutions  to  the  following 
equations  for  wave  potentials  $  and 

I1  '  <V'C/  +  “I  "  0  <A.i> 

35  3x2 

with  a  similar  equation  for  when  is  replaced  by  Cg.  For  any  non-zero, 
constant,  speed  of  propagation,  the  eigen  functions  can  be  derived  to  be: 


U1  "  £  uln0n  5  u2  "  §  u2n0n  ;  °a0  "  S  °agnen  ;  a,&  "  1,2  (A.2,3,4) 


where 


uln  "  u  F(as’ad)  I<n/2)  +  11  frl/2  Cos  (n  °l/2) 
-  05>  g(n)  r®^4  Cos  (n  ©2/2)> 


(A.  5) 


u2n  "  y  F<as’otd)  f<a/2)  +  11  {"“d  rl/2  S±n  (n  0l/2) 

+  (k)  tg(n)/a8]  Sin  (n  ©2/2)>  (A.6) 

°lln  “  F(as»ad)  <n/2)  I(n/2>  +  11  {(2ad  ’  “s  +  X)  r{<n/2)  "  11 
x  Cos  [((n/2)  -  1)  01J  -  g(n)  r*(n/2)  '  1]  Cost((n/2)  -  1)021>  (A.7) 

°22n  “  F(as*ad)  <n/2)  l(n/2)  +  11  {'(1  +  as)  r{(n/2>  “  11 
x  Cos  [((n/2)  -  1)  01]  +  g(n)  rJ(n/2)  "  11  Cost ((n/2)  -  1)  ©2]>  (A. 8) 

°12n  “  F<VV  (n/2>  l(n/2>  +  11  {"2  “d  rl<n/2)  “  11 

x  Sin  [((n/2)  -  1)  0^  +  (h)  1(1  +  ajj/a,!  g(n)  r*(n/2)  *  11 

x  Sin  [((n/2)  -  1)  Oj]  (A.9) 
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V 


where  X,  y  are  Lam&'s  constants;  and  Cs  respectively  the  dilatational  and 
shear  wave  speeds  C.  *  [ (A+2y  )/p]  ;  C  ■  (y/p);  and  the  various  parameters 

Q  S 

appearing  above  are  defined  as: 


a2  -  [1  -  <v/Cd)2]  ;  “3  -  [1  -  <v/Cs)2] 


F<ad’as> 


(1  + 


3(210*  4osad  -  (1  +  a2)2 


(A.  10) 
(A. 11) 


g(n)  *  (4a.a  )/(l  +  cO  when  n  is  odd 

u  S  8 


[1  +  a  ]  when  n  is  even 
s 


10 

tj  e  1  ■  5  +  i  “d*2 


r2  ei01  ■  €  +  iotsx2 


(A.  12) 
(A. 13) 
(A.  14) 


when  v  =  0,  the  above  functions  can  be  reduced  to  the  usual  Williams'  [22] 
eigen  functions. 

It  is  interesting  to  note  that  the  stress  field  Oag(x^)  [oi,B,y  -  1,2], 
should  in  general  case,  satisfy  the  equations: 


n\ 

‘’cB.S  '(  ,tJ  ' 


92u  2  32u 

2v  rrr- —  +  v  - £  p 


3£3t 


(A. 15) 


However,  it  can  be  seen  that  the  special  eigen  functions  given  in  (A. 7-9), 


corresponding  to  the  solution  of  Eq.  (A.l),  satisfy  only  the  equations: 

~2 

,  3  u 

vr»’  i?-° 


(A. 16) 


for  all  values  of  v;  thus,  when  v  ■  0,  the  correspondingly  reduced  eigen 
functions  in  Eq.  (A. 7-9),  which  coincide  with  the  well-known  Williams' 
eigen-function,  needless  to  say,  satisfy  the  static  equations  of  equilibrium, 

°c8,8  "  °* 
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DETAILS  OF  FINITE  ELEMENT  EQUATION  DEVELOPMENT  FOR  DYNAMICALLY  PROPAGATING  CRACKS 
Upon  substitution  of  Eqs.  (1.22  to  44),  into  (1.21),  we  obtain. 


°  -  g  {[(32  kt  +  52  ?T  ‘  §2  +  Si  *T  +  Si  ?T  ’  2l)]  63*2} 

+  ((?  KT  +  ftT  nT  +  rT  mT  nT  4.  ftT  irT 

+  ~2  Js2  +  Z2  5s2  +  ~2  “s2  ~  2s 2  +  5l  ~sl 

+  §1  ill  *  il  ill  -  Sll>  «2 


(B.l) 


where , 


K  -  f  BT  D  B  dv  j  m  -  f  pNTN 
*  u  =  =  =  ~  **v  *  2 

Rn  Rn 

3  2-f  ?T!2dv+/  f?2d* 

VRn  S0n 

2i  ;T  !i d’  +/.  5Tiids 


(B.2,3) 


(B.4) 


(B.5) 


5,2  12  El  O*  +  »  ’ll  ”2  «2>-«  d''  «-6> 

S  vs 

5,1  -I  »i  d’  +  »  *?/„  “2  <»!>•«  d*  -  /  ”2  ?i  ds  (B-7> 


(B.6) 


& 


;,2  ■  5(  »2  52 dv  ••  -.1  -  Hi  Hi 


®.2  *  -2o  v2/v  "2  «2>  •?  dv 

S 


(B.8) 


(B.9) 


?sl  ■  -2»  »l/  “2  <Hl>'S  * 

S 

3,2  •/„  “2  -2  dv  +  ST*  ?2  h  ds 


(B. 10) 


(B.ll) 


3»i  ■/,  ?’  !i  dv  + 


L*  Hi  Hi ds 


fZ'-* 


Now,  the  conditions  of  "least-squares"  matching  of  displacements,  velocities, 
and  accelerations  between  the  singular  element  and  the  surrounding  regular 

0  ** 

elements,  i.e.,  Eqs.  (1.11,12,13)  are  used  to  express  62  and  ^2  in  terms 
of  the  respective  values  q^,  qg2»  qg£  nodes  along  the  boundary  of  Vg.  Thus, 


h  B  hSsi  5  h  m  *i3si  +  Bi  S,i  ;  h  m  hS isi  +  2  SiSsi  +  JiSsi  (B‘12) 

h  =  ~23s2  5  h  "  ^23s2  +  ?23s2  ;  §2  “  *2qs2  +  2  ®23s2  +  ^23s2  (B*13) 

We  note  that  (A^,  and  C^)  and  (A2,  Bg  and  C2)  are  dependent  on  velocities 

of  crack  propagation  v^  and  V2  respectively.  When  Eqs.  (B.12,13)  are  used, 

Eq.  (B.l)  can  be  rewritten  as: 


tv*  in  ip  ip  ip  ip  *p  *¥* 

0  =  n  {<q2  ~  +  ^2  m  “  S2  +  3l  ~  +  3l  ~  ~  Sr  6q2} 
+  lSs2  l?  +  Ss2  "s'  +  3s2  ?sT  -  SsTl  63s2 


(B.14) 


where 


~s  =  [~2  ~s2  ~2  +  ^2  ~s2  ?2  +  ~2  ™s2  ?2] 


(B. 15) 


?s  “  [~2  *s2  ^2  +  2  «2  ~s2  ~2* 


(B.16) 


”s  "  62  ™s2  62 


(B. 17) 


2s  “  tl  <-5.i  5i  -  5.1  §1  "  ;.i  Si  +  ~s2  +  2si> 


(B.18) 


From  Eqs.  (B.14),  Eqs.  (1.45,46)  were  derived.  It  can  now  be  seen  that  both 
.  *  * 

the  singular-element  matrices  K  and  D  are  unsymmetric.  The  "damping"  matrix 

s  s 

*Hp 

o  Is  a  result  of  the  fact  that  the  total  accelerations  of  a  material  point 
s 

in  the  singular  element  depend  on  Bj* 

It  may  be  of  interest  to  note  that  in  the  evaluation  of  Kg2  of  Eq.  (B.6), 
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M  ‘  ' 


the  integrand  will  have  a  singularity  of  the  type  (1/r^)  and  (l/r2> •  Special 
numerical  integration  shcemes  to  evaluate  this  domain  integral  of  Gq.  (B.6) 
directly,  can  be  developed.  Alternatively,  one  can  use  the  observation  that, 
by  definition,  from  Eqs.  (1.21,  B.l) 


?2  =s2  5§2  =£  «\i  +  »<V2  „ 

S  ^ 


32u2  - 

— 2  5  ua)  dv 


(B.19) 


Using  the  divergence  theorem,  Eq.  (B.19)  can  be  rewritten  as: 

32u2 

£  ~s2  %  -f3v  °lj  '-J  8“I  dS  *f„  +  0<V2  S"i  »-20> 


3£ 


the  second  integral  on  the  right  hand  side  of  Eq.  (B.20)  vanishes  due  to  the 
special  property  of  the  eigen  functions  embedded  in  the  singular  element,  as 
explained  in  Eq.  (A. 16).  Thus,  one  can  write  alternatively. 


(B.21) 


wherein,  the  integrand  is  non-singular  along  V  ,  and  no  special  integration 

s 

schemes  are  necessary. 

Likewise,  it  is  seen  that. 


si  six  5§2  +  IXV2  -J  5“I>  -4  “Ij  ^  5“Zi  ds  ‘B-22> 


once  again,  using  the  property,  as  given  in  Eq.  (A. 16),  of  the  eigen  functions 
1 


Kj  in  V^,  and  using  the  divergence  theorem,  we  write 


3*  KT  6g„  -  f  oj,  vj  6u2  ds  -  f  oj  v }  6u2  ds 

~1  =sl  ij  j  i  Jhl  ij  j  i 


(B.23) 


It  can  easily  be  seen  that  the  above  equation  can  be  simplified  to: 


T  T 

«1  E,x  5?2 


5u2  ds 


(B.24) 
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The  above  simplification  is  possible  because:  3Vg  “  pg  +  +  AZ  +  S^, 

where  pg  is  the  interface  of  the  singular  element  with  surrounding  regular 
elements,  and  E^  is  assumed,  without  loss  of  generality,  to  be  free  of  any 
applied  tractions  at  all  times,  and  is  the  ligament  ahead  of  the  crack-tip 
(along  axis)  in  the  singular  element,  where,  for  mode  I  problems,  *  0, 
and  u*  *  0.  The  boundary  integration  as  indicated  by  Eq.  (B.24)  to  evaluate 
K  ,  may  be  more  convenient  that  to  directly  apply  Eqs.  (B.22)  or  (B.23). 
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t  =  2.0  fj,  sec 

re-adjustment  of 
mesh  at  t  =  20  //.sec 

t  =20/isec 


t  =  30  //.sec 


EXAMPLE  :  v  =  1000  m/sec 

At  =  0.2  /isec 
AS=  0.2  mm 


Figure  Captions 


Fig.  1  Schematic  Representation  of  the  Movement  of  the  "Singular-Element" 

Fig.  2  Schematic  Representation  of  Crack  Growth  in  a  Typical  Problem: 

Constant  Crack  Velocity  v  =  1000  m/s;  At  «*  0.2ys;  AI  *  0.2  mm. 

The  mesh  of  regular  elements  around  the  singular  element  is  readjusted 
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